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ABSTRACT 

One  case  of  analytic  data  and  one  case  of  real  data  were  numerically 
integrated  using  a  5-level  baroclinic  primitive  equation  global  model  of 
the  general  circulation.   The  feasibility  of  using  a  mixed  second  order 
and  fourth  order  space  differencing  scheme  to  improve  the  phase  speeds 
of  meteorological  waves  was  examined.   The  results  indicate  that  mixed 
scheme  tends  to  give  a  better  representation  of  the  phase  speeds  than  the 
second  order  scheme  does. 
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I.   INTRODUCTION 

A  baroclinic  primitive  equation  model  using  a  global,  staggered, 
spherical,  sigma  coordinate  system  which  was  originally  written  by 
Dr.  F.  J.  Winninghoff  and  which  is  now  being  developed  by  the  United 
States  Navy  Fleet  Numerical  Weather  Central  (Elias,  1973)  was  used  in 
this  study.   The  purpose  of  this  investigation  was  to  modify  the  hori- 
zontal space  difference  equations  to  a  mixed  second  and  fourth  order 
scheme  with  the  aim  of  improving  the  phase  speeds  of  the  meteorologi- 
cal waves,  as  proposed  by  Williams  (1972).  Real  time  data  (Elias, 
1973)  and  analytic  data  using  an  analytic  spherical  harmonic  stream 
function  (Neatam,  1946)  were  used  as  initial  conditions.   The  use  of 
the  analytic  case  allowed  the  control  of  temperature  and  moisture  dis- 
tributions, predominant  wave  number,  phase  speed  and  wave  amplitude. 
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II.   BAROCLINIC  PRIMITIVE  EQUATION  MODEL 

The  governing  differential  equations,  written  in  vector  form,  are 
similar  to  sets  used  by  Smagorinsky  et  al  (1965),  Arakawa,  et  al  (1969) 
and  Kesel  and  Winninghoff  (1972).   The  integrations  were  carried  out  on 
a  global,  spherical,  staggered,  sigma  coordinate  system  using  the  con- 
servative-type difference  equations  based  on  the  procedure  of  Arakawa 
(1966) .   The  complete  set  of  finite  difference  equations  are  given  in 
Kesel  and  Winninghoff  (1972)  for  the  non-staggered,  polar  stereographic 
grid.   Appendix  A  contains  the  finite  difference  equations  for  the 
spherical,  sigma,  staggered  system  of  equations.  The  heating,  moisture, 
friction  and  diffusion  finite  difference  terms  are  given  in  Kesel  and 
Winninghoff  (1972). 

A.   PRIMITIVE  EQUATIONS 

The  following  are  the  equations  for  the  model  in  spherical,  sigma 
(a)  coordinates. 

The  dimensionless  vertical  coordinate  is: 

a  =  P/tt  . 

The  vertical  velocity  measure  is  defined  as: 

da 
w  =  --=  -a  . 

The  latitudinal  component  of   the  equation  of  motion   is: 


9utt  1  -9(uuTr)        3  (uvttcos8)  , 

3t     ~  "  a  cos  9    l    9X  86  J 


12 


3a  a 


1  ,tt9<J>    ,    RT  3tt,  _      .  _  ,,. 

E   [  -5T  +         -^f]   +  tt  F0  +  D      .  (1) 

a  cos   9        3A  3A  2  u 


The  longitudinal-component  of  the  equation  of  motion  is: 


3vtt  1  , 3  (uvtt)        3(wttcos9)  , 

3t         a  cos  6    l    3A  39  J 


,    ir3(wv)    ,    TTuutan9  ,. 

+    — ^ — -  + TTuf 

3a  a 


1    rir  3$    ,    RT  3tt,  ,„. 

"  7  [     "39  +         39]  +  V  Fy  +  Dv  *  (2) 


The  First  Law  is: 


3TTT  1  r3(TTuT)      .      3CT7VTCQS8), 

3t  a  cos  9    l      3A  38  J 


ir3(wT) 


3a 


RT      ,  ,3tt   ,  1        ru  3tt    .        3tt  cos  9,., 

+  ■„    t-WTT  +   a(-~-  +  q[      "5T  +  v  "5q"  j)J 

c  a  3t       a  cos  9       3A  39 


+  ttH  +  DT   .  (3) 


The  Moisture  Equation  is: 


3qTT  1  , 3(nuq)        3(TTvqcos9)  , 

■3!  a  cos  9    l      3A  39  J 


+     *A<£ll     +  Q7T  +  D      .  (4) 

do  ^  q 
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The  continuity  equation  is: 


3tt  _  1  .3(ut0        SCyTTcosG).        tt  _3w  ,,.,. 

3t       "  a  cos  9    L    3A  89  J  3a  *  W 


The  Hydrostatic  Equation  is: 


3$  =        RT 
3a        "a 


(6) 


In  the  above  equations,  u  is  the  zonal  wind  component;  v,  the 
meridional  wind  component;  tt,  the  terrain  pressure,  T,  the  temperature; 
f,  the  Coriolis  parameter;  $,  the  geopotential;  q,  the  specific  humidity; 
H,  the  diabatic  heating;  F,  the  frictional  stress;  D,  the  lateral  diffu- 
sion of  momentum;  D„,,  the  lateral  diffusion  of  heat;  D  ,  the  lateral 

'  t  q 

diffusion  of  moisture  source/sink  term;  P,  the  pressure;  a,  the  radius 
of  the  earth;  A,  longitude  and  8,  the  latitude. 

Kesel  and  Winninghoff  (1972)  modified  the  pressure  gradient  terms 
as  suggested  by  Kurihara  (1968)  to  reduce  the  inconsistent  truncation 
error  which  arises  as  a  result  of  the  sensitivity  of  the  model  to  varia- 
tions in  the  height  and  character  of  the  terrain.   The  value  of  the  geo- 
potential is  interpolated  at  neighboring  gridpoint  locations  to  the  same 
pressure  surface  as  at  the  computation  point  to  enable  one  to  compute 
to  geopotential  gradients  on  pressure  surfaces  which  are  thus  synthe- 
sized locally  around  each  computational  point.   The  pressure  gradient 
term  in  sigma  coordinates  is: 

RT 

-  V  $  =  -  Va4>  £=■  Vtt  .  (7) 
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Solving  for  the  gradient  of  tt  in  component  form  one  gets, 


3A   RT  l3A  "  9A  J  w 

P    a 


J2L  =  ZL  r 5*  _  ii  i  rQ>, 

30   RT  l36    99  J  v*; 

P     a 


It  is  shown  in  Appendix  A  that  the  right  hand  side  of  equations  (8)  and 
(9)  are  the  sum  of  the  thicknesses  between  the  sigma  surface  and  the 
pressure  surface  to  the  right  and  left  of  the  computational  points. 
Since  the  thickness  depends  only  upon  the  mean  temperature  of  the  layer 
in  question,  the  thickness  equation  is  used  to  solve  equations  (8)  and 
(9) .   The  mean  temperature  is  obtained  by  interpolation  as  indicated  in 
Appendix  A. 

B.  GRID 

The  longitudinal  and  latitudinal  grid  increments  were  ten  degrees 
for  the  mass  and  velocity  fields.   The  geopotential  ($),  temperature 
(T) ,  specific  humidity  (q) ,  and  the  vertical-velocity  (w  =  -a)  were 
carried  at  the  poles.   When  $,  T,  q,  w  and  it  were  needed  at  a  velocity 
point  they  were  taken  as  the  average  of  the  four  surrounding  mass  points. 
Figure  1  represents  the  staggered  grid  and  location  of  the  variables. 

C.  VERTICAL  LAYERING 

Figure  2  is  the  vertical  layering  as  described  by  Elias  (1973). 
The  basic  variables,  as  indicated,  are  carried  at  the  center  of  each 
layer  except  that  the  specific  humidity  (q)  is  carried  at  the  lower  three 
layers  only.   The  vertical  velocity  (w)  is  carried  at  the  interface 
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Figure  1.   Location  of  Variables  and  Grid. 
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Figure  2.   Vertical  Layering 
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between  layers.   The  vertical  coordinate,  sigma  (a)  was  defined  by 
Phillips  (1957). 

D.   TIME  DIFFERENCING 

The  two.  finite  difference  techniques  that  were  used  were  the  cen- 
tered difference  scheme  and  the  Matsuno  (Euler-backward)  difference 
scheme.   The  finite  difference  equation  for  the  centered  technique  is 


Ft+1  =  F^1  +  2At  |£ 

ot 


and  that  for  the  Matsuno  scheme  is 

Fk   =  F*  +  At  9Ft 


9t 


Ft+i  =  Ft  +  At  |£ 

ot 


Since  the  centered  scheme  is  not  feasible  for  the  first  time  step  and 
also  produces  a  computational  mode  (Haltiner,  1971),  the  Matsuno  scheme 
was  used  for  the  initial  time  step  of  each  3-hour  integration  to  reduce 
solution  separation.   The  Matsuno  scheme  also  selectively  dampens  high 
frequency  waves  (Haltiner,  1971).   The  Arakawa  technique  (Langlois  and 
Kwok,  1969)  of  averaging  quantities  involved  in  the  longitudinal  deri- 
vitives  is  used.   This  allowed  a  time  step  of  ten  minutes  to  be  used. 
Without  this  technique,  the  Von  Neuman  linear  computational  stability 
criterion  (Haltiner,  1971)  would  require  a  two  and  one  half  minute  time 
step  because  of  the  longitudinal  grid  distance  at  85°  latitude  which  is 
only  47  kilometers. 
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The  lateral  diffusion,  friction,  convective  condensation  and  large 
scale  condensation  are  computed  at  every  time  step  while  heating  is  com- 
puted every  six  time  steps  (Kesel  and  Winninghoff,  1972). 

E.   FOURTH  ORDER  MIXED  SPACE  DIFFERENCING 

Williams  (1972)  found  that  the  stability  criterion  for  fourth  order 
unstaggered  space  differencing  required  a  time  step  which  was  27%  smaller 
than  that  for  second  order  space  differencing.   He  then  noted  that  the 
phase  speed  of  the  meteorological  wave  is  independent  of  the  difference 
approximation  which  is  used  for  the  pressure  force  term.  A  mixture  of 
fourth  order  differencing  in  the  advection  terms  and  the  divergence  term 
in  the  continuity  equation  and  second  order  differencing  in  the  pressure 
force  term  was  proposed.   The  stability  criterion  for  this  mixed  scheme 
was  found  to  require  a  time  step  that  was  only  14%  smaller  than  that  of 
the  second  order  scheme.   Thus  the  mixed  differencing  allows  a  larger 
time  step  while  giving  fourth  order  accuracy  in  the  phase  speed  of  the 
meteorological  wave.   Appendix  B  presents  the  fourth  order  differencing 
of  the  advection  terms  and  the  divergence  term  of  the  continuity  equa- 
tion.  A  time  step  of  8  minutes  was  used  for  this  scheme. 
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III.   INITIAL  CONDITIONS 

Two  cases,  one  using  real  data  and  the  other  analytic  data,  were 
used  for  initialization  of  both  the  second  order  and  mixed  differencing 
schemes.   The  real  data  was  produced  external  to  the  main  program  on  a 
Northern  Hemispheric  FNWC  63  x  63  grid.   The  analytic  data  with  known 
phase  propagation  was  produced  following  the  work  of  Neamton  (1946), 
Gates  (1962),  Heburn  (1972)  and  Elias  (1973).   The  real  data  were  fields 
analyzed  by  FNWC  objective  schemes. 

The  main  program  was  used  to  interpolate  the  data  so  that  values  as 
5°  latitude-longitude  intersections  were  available  on  the  Northern  Hemi- 
sphere (Elias,  1973).   The  data  was  then  reflected  into  the  Southern 
Hemisphere.   In  both  cases,  the  input  data  consisted  of  temperature 
analyses  for  the  Northern  Hemisphere  at  12  constant  pressure  levels  from 
1000  to  50  mb,  height  analyses  at  ten  of  these  levels,  moisture  analyses 
at  four  levels  from  the  surface  to  500  mb,  sea  level  pressure  and  sea 
surface  temperature.   Monthly  mean  surface  temperature  fields  were  used 
to  derive  an  albedo  field  (Dickson  and  Posey,  1967).   In  both  cases  the 
terrain  height  was  set  to  zero.   For  the  barotropic  analytic  case  heat- 
ing, convective  adjustment  and  friction  were  turned  off. 

If  one  assumes  nondivergent  horizontal  flow,  harmonic  wave  solutions 
of  the  complete  vorticity  equation  can  be  obtained  for  the  sphere 
(Neamton,  1946).   These  solutions  give  the  velocity  of  propagation  of 
the  wave.   The  solution  of  the  \p   field  (Elias,  1973)  was  found  to  be 

\b   =  A  sin  (mA-vt)Pm  (sin  9)  -  B  a  sin  6  +  C  P   (sin  9)      (1) 

n  n 


20 


where  A,  B  and  C  are  constants  to  be  determined,  a  is  the  radius  of  the 

earth,  m  is  the  hemispheric  wave  number,  v/m  is  the  angular  phase  speed 

of  the  wave  (radians/second),  P  denotes  the  Legendre  polynomial  and  P 

n  n 

represents  the  associated  Legendre  function. 

If  C  =  0  and  n  =  m  +  1  (Elias,  1973),  equation  (1)  reduces  to 

\\i   =  A  sin  (mX  -  vt)  (2N!/2NN!  )sin  9  (cos  9)6  -  B  a2  sin  6      (2) 

The  constant  A  is  arbitrary  and  proportional  to  the  wave  amplitude.   It 
was  shown  by  Haurwitz  (1946)  that  the  solution  obtained  for  ^  implies 
the  existence  of  a  velocity  distribution  over  a  sphere  such  that  the 
angular  velocity  of  the  westerly  current  is  made  up  of  the  sum  of  three 
terms.   The  first,  B,  is  constant  over  the  sphere;  the  second  term 
varies  along  a  meridian  but  is  constant  around  a  latitude  circle;  and 
the  third  term  represents  the  angular  phase  speed  which  is  given  by 

V/m  =  B  «¥&£   -  §°  (3) 

n(n+l)    n(n+l) 

where  ft  is  the  earth's  rotation.   For  the  analytic  experiment  A  was  chosen 

2   -1  -7   -1 

as  1000  m  sec  ,  m  as  6,  n  as  7  and  B  as  8  x  10  sec   .   The  solution  of 

equation  (  3)  gives  an  angular  velocity  of  -9°  longitude/day  (toward  the 

west).   Since  wave  number  six  was  used,  integrations  were  performed  for 

only  one-sixth  of  a  latitude  circle  with  the  remaining  five-sixths  being 

set  equal  to  the  first  part.   This  resulted  in  a  40%  saving  of  computer 

time  for  the  analytic  case  only. 

The  geopotential  fields  were  derived  by  solution  of  the  linear 

balance  equation  and  the  temperatures,  constant  at  each  pressure  surface, 
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were  consistent  with  the  National  Advisory  Committee  for  Aeronautics 
(NACA) .   Both  the  geopotential  and  temperature  fields  were  derived  as 
shown  by  Elias  (1973). 

The  wind  field  was  produced  following  the  work  of  Elias  (1973)  by 
using  the  finite  difference  expressions 

u  -  ~  ~     A9  <*> 


*  (5) 


a  cos  9  AX 

where  a  is  the  earth's  radius,  A9  is  the  latitudinal  distance  increment 
and  AX  is  the  longitudinal  distance  increment.   Poleward  of  25°  latitude 
the  initial  geopotential  field,  $  ,  was  retained  and  the  ij;-field  was 
obtained  by  the  solution  of  the  linear  balance  equation, 

2  2 

V  4>  +  Vip   •    Vf/f  »  V  <J>   /f  (6) 

Equatorward  of  25°   latitude  the  infield  was  determined  using, 

$  -  *  /! 

where  f  is  a  mean  Coriolis  parameter.  A  new  geopotential  field  ($) 
was  obtained  using 

2         2 
V  *  =  (fV  I})  +  Vl^r  •  Vf)  . 

At  25°  latitude,  ifi  and  ty,   were  a  combination  of  half  the  value  of  the 
poleward  case  and  half  of  the  equatorward  case. 
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IV.   WAVE  ANALYSIS  METHOD 

A  Fourier  series  was  determined  at  each  five  degrees  of  latitude 
around  the  latitude  circle.   With  this  technique,  the  phase  angles  and 
amplitudes  of  each  wave  number  around  a  latitude  circle  can  be  calcu- 
lated. A  Fourier  series  can  be  expressed  as  follows  (Heburn,  1972): 


F(x)  =  A  +  >^  (A  cos  mx  +  B  sin  mx) 
o   Z^    m  m 

m 

=  C  +/  (C  cos  (mx  -  6m)) 
o  *—*  m 


m 


B         A 
where        C  = 


"m   sin(6m)    cos(6m) 


and 


6m  =  tan-1  (W 

m  mm 


The  analytic  experiment  involved  an  input  stream  function  of  wave 
number  six.   The  values  of  primary  interest  were  6,  and  C,  which  are 
the  phase  angle  and  amplitude  of  wave  number  six.   Other  angles  and 
amplitudes  were  extracted,  especially  in  the  real  data  case. 


23 


V.   RESULTS 

Two  experiments  will  be  presented  in  this  section.   The  first 
experiment  was  performed  with  analytic  fields  and  the  second  with  FNWC 
analyzed  fields  from  1200Z  10  May  1973. 

As  indicated  in  section  III  the  angular  phase  speed  for  wave  number 
six  was  -9°  longitude/day  (toward  the  west).  A  72-hour  forecast  of 
surface  pressure  was  made  by  both  the  second  and  mixed  order  differenc- 
ing.  The  changes  in  phase  angle  vs  latitude  were  plotted  for  both 

i 
cases  and  are  shown  in  Figures  3  and  4.   Table  I  summarizes  the  average 

angular  movement  per  24-hour  the  entire  72-hour  forecast  for  both  dif- 
ferencing schemes.   It  should  be  noted  that  the  mixed  differencing 
scheme  propagated  the  wave  slower  to  the  west  than  the  second  order 
scheme.   The  barotropic  Rossby  equation  (Holton,  1972)  is  given  as 

C  =  u-    * 


2   2 
(k  +  hi  ) 


where  3  is  -j—  ,  f  is  the  Coriolis  parameter,  u  is  the  mean  zonal  flow, 
k  is  the  zonal  wave  length,  m  is  the  wave  number  and  C  is  the  zonal 
phase  speed.   Since  the  quantities  3>  k  and  m  are  constant  the  phase 
speed,  C,  depends  only  upon  the  mean  zonal  flow.   Fourth-order  differ- 
encing should  give  a  better  approximation  for  the  effect  of  the  mean 
zonal  wind,  u,  and  therefore  the  westward  movement  of  the  wave  would 
be  slower  than  that  of  the  second  order  differencing  case.   It  would 
appear  that  the  mixed  differencing  scheme  does  give  a  better  solution. 
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Figure  3.   Phase  Angle  vs  Latitude  for  the  Second  Order  Scheme. 
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Longitude 


Figure  4.   Phase  Angle  vs  Latitude  for  the  Mixed  Scheme. 
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SECOND  ORDER 

MIXED  SCHEME 

LATITUDE 

DEG  LAT/DAY 

DEG  LAT/DAY 

70° 

-12.6° 

-4.8 

60° 

-10.3° 

-7.9 

50° 

-9.2° 

-7.4 

40° 

-8.5° 

-7.5 

30° 

-7.7° 

-7.2 

20° 

-7.3° 

-7.6 

10° 

-7.0° 

-6.5 

0° 

-4.6° 

-3.6 

TABLE  I.   Average  Movement  of  Analytic  Wave  per  24-Hours  for  a 

72-Hour  Forecast  for  Both  Second  Order  Differencing  and 
the  Mixed  Scheme. 
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A  48-hour  forecast  was  also  made  by  both  the  second  order  and 
mixed  differencing  schemes  using  analyzed  FNWC  data  fields.   Chart  A 
is  the  initial  analyzed  surface  pressure  field,  Charts  B  through  E 
are  the  12-,  24- ,  36-  and  48-hour  forecasts  produced  by  the  second 
order  scheme  and  Charts  F  through  I  are  the  corresponding  forecasts 
produced  by  the  mixed  differencing  scheme.   Table  II  lists  five  high 
pressure  centers  and  Table  III  lists  seven  low  pressure  centers  with 
their  initial  position  and  both  the  48-hour  second  order  and  48-hour 
mixed  scheme  forecast  positions.  As  can  be  seen  each  high  pressure 
center  was  forecast  at  least  five  degrees  farther  to  the  east  by  the 
mixed  differencing  scheme.   Low  pressure  centers  1  and  3  also  were 
forecast  five  degrees  farther  to  the  east  by  the  mixed  scheme  while 
low  centers  2  and  4  were  moved  very  little  by  both  differencing  schemes. 
Low  number  5  was  the  only  low  that  was  moved  slower  than  the  second 
order  scheme.   Low  number  6  was  developed  into  two  centers  by  both 
forecast  schemes.   Only  the  forecast  position  of  the  major  center  was 
listed  in  Table  III.   The  mixed  scheme  moved  this  center  nine  degrees 
farther  than  the  second  order  scheme.   Center  number  7  was  also 
initially  a  single  center  but  was  developed  into  three  small  centers 
with  the  second  order  forecast  scheme.   The  mixed  scheme  maintained 
the  more  realistic  single  center  system  near  its  initial  position. 
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Chart  A.   Initial  Surface  Pressure  Field. 
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Chart  B.   12-Hour  Forecast,  Second  Order  Scheme. 
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Chart   C.      24-Hour  Forecast,    Second  Order  Scheme. 
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Chart  D.   36-Hour  Forecast,  Second  Order  Scheme. 
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Chart  E.   48-Hour  Forecast,  Second  Order  Scheme. 
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Chart  F.   12-Hour  Forecast,  Mixed  Scheme. 
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Chart  G.   24-Hour  Forecast,  Mixed  Scheme. 
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Chart  H.   36-Hour  Forecast,  Mixed  Scheme. 
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Chart  I.   48-Hour  Forecast,  Mixed  Scheme. 
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HIGH  CENTER 

INITIAL 

48-HOUR  SECOND 

48-HOUR  MIXED 

NUMBER 

POSITION 

ORDER  POSITION 

SCHEME  POSITION 

•  1 

37N  177W 

37N  17 7W 

35N  172W 

2 

35N  13 6W 

47N  128W 

43N  121W 

3 

None 

31N  53W 

33N  48W 

4 

45N  25W 

52N  17W 

52N  12W 

5 

64N  50E 

60N  60E 

52N  66W 

TABLE  II.   Initial  Positions  of  High  Centers  and  Forecast  Positions 
of  the  Second  Order  Scheme  and  the  Mixed  Scheme. 


,0W  CENTER 

INITIAL 

48-HOUR  SECOND 

48-HOUR  MIXED 

NUMBER 

POSITION 

ORDER 

POSITION 

SCHEME 

;  POSITION 

1 

47N 

155W 

48N 

152W 

49N 

14  7W 

2 

47N 

89W 

50N 

45W 

48N 

40W 

3 

47N 

50W 

50N 

87W 

51N 

88W 

4 

67N 

02E 

68N 

03E 

68N 

04E 

5 

6  IN 

89E 

59N 

99E 

61N 

97E 

6 

36N 

152E 

50N 

151E 

48N 

160E 

7 

30N 

80E 

* 

33N 

81E 

TABLE  III.   Initial  Positions  of  Low  Centers  and  Forecast  Positions 
of  the  Second  Order  Scheme  and  the  Mixed  Scheme. 

*Several  small  lows  were  forecast  by  the  second  order  scheme. 
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VI.   CONCLUSIONS 

Both  the  analytic  and  real  data  experiments  indicate  that  the 
mixed  fourth  and  second  order  differencing  scheme  gives  a  more 
accurate  solution  for  the  meteorological  phase  speeds.  More  analytic 
experiments  with  a  more  realistic  phase  speed  will  have  to  be  per- 
formed to  produce  the  final  conclusive  proof  of  the  mixed  scheme. 
A  series  of  real  data  experiments  should  also  be  considered.   In 
addition  the  fourth  order  differencing  must  eventually  be  extended 
to  the  poles. 
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APPENDIX  A 
DIFFERENCE  EQUATIONS 

The  longitudinal  momentum  equation  terms  in  difference  form  are: 
n+1     n-1 

iHl  =    1-1         XJ  (1) 

9t  -      2At       '  U; 


1  , 9  (uutt)        3(uvtt)cos81 

a  cos  9    l   3A  96  J 


n  (U-.+U...      .)u7T.,1/0      .     -     (U..+    U.     ,      .)uTT.     ..    ,- 

1        f    xj       i+l,3       i+l/2,,i  ij       i-i, y      i-i/2t.i 

a  cos  6. l  2AA. 


(u..+  u.     .,0)vtt.     ...cos   9  -    (u..-  u.    .    0)vtt.     .      COS   6.    -     (2) 

+        13        i,j+2        i,j+l  j+1  ij        i,j-2        i,j-l  j-1 


where 


Uir.L1/,      .     =    7(l..U..+    I.j,      -U.,1      .+    IT-      .,,U.      .,-,+    TT.       ,     ,U.      ,     n)      , 

1+1/2,3        4      ij    ij        l+l, j    l+l, j        i,J+l   i,J+l        i,3-l   i,J-l 


(2.1) 


uiT.    n  ,_    .   = -r(i.  .u. .+  ir.   n    .u.    n    .+  tt.    .    ,,.u.   ..    ...    . 

i-l/2, j       4V  ij   13       i-i, j   i-l, j       i-i, j+1  i-l.j+1  + 


ir.    ,    .    ,u.    ..    .    ,)  (2.2) 

l-l.J-l   1-1,3-1 


UTT.      ...  =   7-(lT.  .U.  .+    TT.      .,„U.  +   7T  >iu-     t      -,i    + 

1,3+1  4      13    ij         1,3+2    1,3+2        1-1,3+1   i-l,J+l 


'l.J-A.J-l*  (2'3) 
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U^4   4   1   =  t(™1U-+    ^4   •   0U4   •   0+  ^4   1   4   1U4   1   4   1   + 

i,j-l   A  ij  ij   i,J-2  i,j-2   i-l,j-l  i-l,j-l 


'iJ-Aj-i*  (2'4) 


The  longitudinal  grid  increment  (AX)  and  the  latitudinal  grid  increment 
(A0)  are  both  ten  degrees. 


»         w.  .,  (u.  .  ,  , .,  +  u.  .,  )   w.  .  .  ,  (u.  .,  +  u.  .  ,  n) 
tt  9wu  =  -   ,  ilk  u.k+l    ink  _  i,i,k-l  i.ik n,k-l     m 

11  3a  "  ijL      2Ao  2Aa  ^J; 


where  the  vertical  velocity,  -a  ,  is  defined  as, 


Aq        f3ir>H 
Wijk  "  -        |_C8t;ij 


Wk+1     -  "ijk      -      rat'ij 


.  UTT..T      .,      -    UlT  .,      -    UTT...  V7T.  .  COS0  .  .  , -VTT  .  .     ..COS.. 

+  1  r      i+l>.lk  i+l,  ,ik  i.i  k  ii+lk  .i+l        i.i -Ik  .1- 

a  cos   6.  L  AX  AG 

J  (4) 


"] 


and   the   local   change  as 


^       «  5        UlT  -UTT...  VTT.  .,-.  COs8.  ,-,-VTT.  .     nlCOs8.     , 

,3n\N  _  V"^  r     1+1.1  k       i.i  k  13+lk         .i+l       i.i -Ik         j-1. 

Wi-rZ-f      AXa  cos   9.  A6  a  cos   6.  J  w 


k=l  J  d 


The  utt   terms   of   Equations    (4)   and    (5)   are, 


UTT.  .,  .     =    U.  .,  .  (TT.  .+   TT.  ,  ,  .+   7T.  .  .,+   IT.  .     ,)   t  C5.1) 

i+l  j  i+lj      ij        i+lj        ij+1       ij-1     4  w*    ' 


The  quantities  urr .  . ,    utt..,.    and  utt  .  .    ,    are  also  determined   in   the 
ij'        ij+1  ij-1 

same  manner. 
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Tfvf  =  f  ^i^H  (6) 


where 


ij   4V  ij   l-lj   l-l, J +1   i-l, J -1' 


_   u .  .  (irv .  . )  tan9 . 
Tfuv  tan8    ij    ij J_ 


1    rRT  3ir  .  it  _8J>, 


(7) 


a  cos  6     9X     9X 


1 


[RT, ,  fe*  +  i,  -JV-fcL*-.]  (8) 


a  cos  6     ij  9X      ij    AX 


where 


f..  =  y(T.  .,+  T.  .+  T.  !  .,,+  T.  .  .  ^)  (8.1) 

ij   4V  l-lj    ij    i-l, j+1   i-l, j -1 


ir. .  =  tO*.  ,,+  ir..  +ir.  ,  ...  +  ir.  .  ,  .,)  (8.2) 

ij   4V  i-Ij   ij    i-l, J+1    i-l, J -1 


From  Equation  (9)  of  section  II  the  change  of  terrain  height  is  found 
as  follows 


,9*   *    7T    <»R-  V  <»R  -  V.  fQ, 

(3Aa>  =  if  C   AX    P AT—  a]  (9) 

<*R  -  *R  )+  <*L  "  *L  > 

=  Z_  r E 2 2 E_i  (10) 

RT  l  AX  J  K±UJ 
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where  $  is  to  the  right  of  the  computational  point  and  <J>  is  to  the 
left  of  the  computational  point.   Sigma  and  p  indicate  the  surface  in 
question.   The  thickness  between  the  pressure  surface  and  the  sigma 
surface  to  the  right  and  left  of  the  computational  point  is  defined  as, 


MR   =  <$R  "  $R  }      and  (11) 

Op       p      0 


A$l   =  ($L  -  $L  )  (12) 

ap     a    p 

therefore  one  can  write, 


op      "r    v2 
The  mean  temperature  is  found  as  follows, 


Pi 

A<J>n       =  -RT       In  —  (13) 


Pi 
A$T        =  -RT        £n  —  (14) 

L^  Op  p. 


comp       da       2 


and   in  difference   form 

(TR   ,      "  TR  '  (LN  TT .  .    -  LN  tt  .  .  ) 

h      -  hi  +        °lo        °      S *""  <16> 

ap  J 

(TL    . .    "  TL   )  (LN  7T .  .    -  LN  TT .  .  ) 

f          =  f        +  5+i 1 U U_  (17) 

lL            Xij                   Aa  2                                                         u/; 
aP          J 


A3 


Vj  =  i(Tij  +  Ti-ij  +  Ti-i,j+i  +  'i-ij-i*  (17-1) 


LN  if .  .   =  rC^n  ti\  .   +  Jin  tt.    ..  .   +  Jin  tt.    n     . .,  +  An  IT.   ,     .    n)  (17.2) 

ij        4V  ij  i-lj  i-l,j+l  i-l, j-1' 


Figure  5  illustrates  the  differencing  in  both  the  horizontal  and  verti- 
cal. Substituting  Equations  (16)  and  (17)  into  Equations  (13)  and  (14) 
one  obtains, 

(TR   .,"TR       (Jin  tt.  .-£n  tt..) 
A\p  "  -RtV  AC       ° 2 ^^  \fZn  V      (18) 

(TL         ~TL    )        (£n    ^_-ln    w  ) 

A$T        =-R[T..4    — 2±i 2_    11^ ^J_]    £,„tt .  . -Jlmr .    ,.      (19) 

Lap  l   ij  Aa  2  ij  i-lj 


where 


P-i  (a+1    Att) 

j>n     -^ —  =   Jin  — — =   Jin  TT.  .    -   Jin  TT.  .  (20) 

^Pa+1  W  1J  1J 


and 


Pi 
JinT   -= —  =   Jin  tt.  .   -  Jin  tt.    .  .  (21) 

L  PG+1  1J  X-lj 
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Figure  5.      Interpolation  for  Vtt. 
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The  latitudinal-momentum  equation  terms  in  difference  form  are, 

N+l       N-l 

•>      TV.  .      -  TTV.  , 

Ivrr  =   ij ij  (1) 


3t  "       2At 


1  r  3  (ttvw)   3(Trwcos8)  ,    _ 


a  cos  8  '   8X        36 


-  (V..+V  )UTT  .  -     (v..+V.         .)uTT.         . 

1  r         IJ    •     1+lj  1+1/2,1 1,1         1-1,1  1-1/2J 

a  cos  6.  L     ~*~  2AA 

J 


(V.  .+V.  .  ,0)VTT.  .  ...COS6.  -     (v.. -V..     0)VTT..     ^OsB.     .. 

.        11      iJ+2        i.i-fl  1+1  1,1      1.1-2        i.i-l jjvL,  m 

+  2A6  J     K  } 


The  utt  terms  are  determined  in  the  same  manner  as  the  terms  for 
Equation  (2)  were  found. 


r.  /      \  w    (v .  .    , ,    +  v .  .    )  w      ,  (v .  .      +  v .  .      .) 

tt3(wv)         -     ,  av  i,ia+l         ijcT  a-lv  13  a         13  a- 1  \  ,,, 

3a    "  ij L     2Aa        "       2Aa       J     u; 


The  vertical  velocity,  w,  and  the  local  change,  (-5—)  .  .  ,  are  obtained 

dt  xj 

as  in  Equations  (4)  and  (5)  of  the  longitudinal  equation  of  momentum. 


TTvf  =  f .  tt.  .  v.  .  (4) 

J   13   13 


IT..  =  Equation  (6.1)  or  the  u-equation  of  momentum  + 


-      .   u.  .  (tt.  .u.  .)tan6 . 

tuu  tan9  _   H   H  H      1  ( 5) 
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1   .it  3<t>       RT  d-n, 
"  a   l     39  36J 


1      ds„ (U>*  +  i„  (*1-1"1*1,:  *l-l-1'l)]  (6) 


a  cos  9      l     ij    39;,  ij  A9 


<*R     -  *R   >   +   <*L     "   *L   > 


<W>*     '  S  f-« V"^ «-]  .  <7) 


Equation  (7)  is  solved  by  the  same  method  as  was  shown  for  the 

lontitudinal  momentum  equation. 

Thermodynamic  Equation 

N+l    „   N-l 
-  ,  mN   ttT  .  .     -  ttT  .  . 

3t   "        2At  v  ; 


1     r3(TruT)     3(TrvTcos9  . 


a  cos  9   L   3A         39 


.  T...  .(u.  .,  .It...  .)  -  T.  .(u.  .TT.  .) 

1  r  i+lj   1+1,1  i+l J      ij   ij  ij 

a  cos  9.  AA 
J 


T...-(v.  ...IT..  .1)cos9.,-  +T..  .v..  ,17..  _cos9.  . 
x,n+lx  ij+i  i,i+l    jjjL    ij-i  l.i-i  ij-i   j-lj  (2) 


where  T  is  the  temperature  at  the  velocity  point  based  on  the  average 
of  the  four  surrounding  mass  points  and  tt  is  also  calculated  in  the 
same  manner. 
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*  Sfr?)       .       rwi1o(Ti,ia+l  +  Tijo)  wiia-l(uija  +  "ijo-l*         m 

9a      ~  uij    l  2Aa  2Aa  J        u; 


The  vertical  velocity,   w,   and  the   local  change,    (-^— )      ,    are  obtained 


as   in  the  u-equation  of  motion. 


RT 
c 


[,  ,377  .   1    ,  377  ,   377    alNl 
-wTT+a(-vr  + Q[u  "tt  +  v  ^q  cos9J) 
3t   a  cos  8    3A     86        I 


RT  r(wq  +  wq-l)  w     n     ,9tt. 
c^  [ 2 *U  +  °a  <5t>iJ  + 


a     (u.11.ir.J1.-u..Tr..) 


a  cos  6.  l         2AA 


(v. .  ...77. ..  ncos8.  ,  n    -   v..  .TT..  ..cosB.  ..) 

xj+1  ij+1    j+1    ij-1  i.T-1 .1-1  ,  ,,, 

2A6  J  ^ 


where  the  vertical  velocity  and  the  local  change  are  again  computed 
as  before. 


The  Moisture  Equation 

.     N+l  N-l 

qtt .  .  -  qTT .  . 

(1) 


.     N+l  N-l 

9  qTT    _        ij ij 

3t       ~  2At 


1  ,duTrq    ,    3(v77qcos8) .      _ 

a  cos  8      l    3A  38  J      ~ 

,  q.  ,,  .  (u.  ..  .TT.  ,  ,  .)  -  q..(u..77..) 

I r  1+1-1   1+11  1+1-1 1.1   ij  U      + 

a  cos  8  .  AX 

J 

q.  .  .,  (v.  .  ,,tt.  .  .  n)cos6  .  ,  -  -  q..  ,  (v  .  ,77..  -cos6.  .. ) 

Hx2+iy  i,i+i  jj+r       i+i     Vi-iv  n-i  ij-i       ,-i-r, 

Aa  J 


(2) 
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where  q  is  the  specific  humidity  at  the  velocity  point  based  on  the 
average  of  the  four  surrounding  mass  points  and  tt  is  also  calculated 
in  the  same  manner. 


it  3(wg)      rwija(qiio+l   qijcr 
do  ij l       2Aa 


w.     (u.   +  u.._  ,) 

2Aa       J  VJ; 


The  Hydrostatic  Equation 

d<\>  =  -  RT  d£n  O 


\  ak+l 

ak+i  Vn(-¥>  +  Wn(  *    > 

$,,.,=$,     -   R£n  [ 5 ; — : — r^ ] 

k+1         k  a,  £na     +  £na         -  2£niT   J 


At  the  lowest  level  $  is  given  as 


*..  =  TT  -  RT.  ..  In    (0.9) 
ij         ijk 
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APPENDIX  B 

In  accordance  with  the  work  of  Williams  (1972)  the  fourth  order 
terms  are  a  sum  of  a  weighted  average  of  the  second  order  and  fourth 
order  differencing.   The  weighting  factor  the  second  order  part  is 
four-thirds  and  a  minus  one-third  for  the  fourth  order  portion.   The 
advection  term  of  the  longitudinal  momentum  equation  is: 

.  .         (u.  .+    U.  ,  ..  .)uTT.     1/0.     -     (u..+    U.     ,  .)uTT.     1/0. 

1  r4      rV     1,1  1+lj  1-1/2J  1J  1-lj  1-1/2J 

"  a  cos  6.    l3  2AA 

J 

(u.  .+  u.  .  .  _)cos8_  1   VTT,..,    -(u..-u..    „)cos9.    n    vrr.. 

2AA  } 


n    (u.  .+  u_  0.)uir.      .-(u.  .+  u.    „.)utt. 
_  V   i.l        i+2,l        i+l,l        Li        x-2.1        l-l.l 

3l  4AA 


(u .  .+  u.       . )  cos8      -   VTT  .  -  (u  .  .+  u .  .    .  )  cos9  .         VTT  .  . 

AAA  n      U; 


The  viT  terms  are  determined  as  shown  in  Appendix  A.   The  advection  term 

for  the  latitudinal  momentum  equation  is  similar  to  term  (1) .   The 

advection  term  of  the  thermodynamic  equation  is: 

n  .  f ....  .    UTT.  ..  .  -  T.  .    uif .  . 

_         I  A  r   1+1-i       1+1-i         Li LL 

a  cos  6.  l3  AA 
J 

T.  .  ^utt.  _,cos9.  .-    -  T..    .    uif..    .cosG.    n 
i.l+l     l.i+l         i+l         i.i-l       l.i-l 2zL\ 

AA  ' 

_*         ~^*  -*         —=* 

1     T.,..     UTT  .     -    T.         .     UTT.     . 

_  If   i+l.l        i+1.1  1-1.1 izil  + 

3X  2AA 

— *     — *  _£     ~* 

I...,  UTT.    cos6.  ,„  -  T..  _  UTT..  0cos6.  _ 
ij+2   ij+2 j+2 ij-2 ij-2 1^2,,  m 

2A9  ;J  v  ; 
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The  program  as  written  gives  the  average  temperatures,  T,  and  utt  values 

at  velocity  points  for  the  second  order  differencing.   For  fourth  order 

_*      _* 

differencing  the  grid  distance  is  doubled  and  the  values  T  and  utt  are 

determined  at  the  adjacent  mass  points.   All  the  averaged  quantities 
are  determined  similarly  to  those  in  Appendix  A. 
The  flux  term  of  the  continuity  equation  is: 

.   utt.  A.  .,    -   utt.  ., 


a  cos  8.  L3       AX 


vir. . , ,,  cos0. , .  -  vrr..  n,cos8.  , 
n+lk    j+1     it -Ik    .tIi 

AX  } 


~*       — —* 
3  l      2AX 


V7T.  .  -,  COS0.  -  -  VTT..  „COs9.  „ 

ij+2k    j+2     n-2    ]-2n 
2A9  n 


The  averaged  values  are  also  determined  as  indicated  in  Appendix  A. 
The  unstarred  quantities  are  at  a  velocity  point  and  the  starred 
quantities  are  at  a  mass  point.   For  this  experiment  the  fourth  order 
differencing  was  blended  with  the  second  order  differencing  at  75° 
latitude.   This  was  done  to  avoid  having  to  go  over  the  poles  to  pick 
up  the  j+1,  j+2  and  j+3  quantities  at  the  higher  latitudes. 
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